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ffiP  qt'iie.  ?il  P';  tri+ -fiiit  'leiCribe  llie  'low  c'  the  ocean  ore  well 
tnowii,  tiul  for  Nan.  ca^ei  of  interest,  analytical  solutions  are  not 
possible  because  of  irre-iular  boundaries  or  significant  no n - 1 ineor it ies 
in  tbe  equations.  It  is  often  possible  to  obtai/i  useful  aporoxiNute 
solutions  usinq  lumerical  techniques.  However,  for  reasonably  detailed 
nuNerical  Nodels.  very  larqe  anounts  of  coMputci  tine  are  required, 
and  it  becoMes  iMperative  to  seek  out  the  nost  efficient  nunerical 
alqoi'ithns.  In  the  case  where  forecasts  are  beinq  Made  on  a  real  tine 
basis,  the  use  of  non-of'linun  schenes  nav  result  in  unacceptable 
dr  lavs . 


fbis  report  descrites  a  class  of  techniques  for  the  efficient  treatment 
of  podies  of  water  witti  irregular  boundaries.  Realistic  nodels  of  the 
ocean  state  nust  include  adequate  treatnent  of  the  boundaries,  the 
irreqular  bcttoM  topoqraphy,  the  shape  of  the  coastline  and  the 
co-Npl  1  cat  pd  qeciNPtrv  of  qroups  of  islands. 

Traditional  nuMerical  Methods  usinq  a  rectanqular  finite  diffc-rence  qeid 
are  suitable  for  exploratory  surveys  of  idealised  pioblens  where  Itie  aiM 
is  to  see  why  a  particular  svsteM  responds  as  it  does,  and  wl.nl  effects 
var ;  at ;  u  ."'S  in  the  pai'aMeters  of  the  p'obleM  and  the  fu.'cinq  terms  tiSive 
on  tue  outcoMe,^  In  order  to  use  these  Models  for  compile  ate.j 
q^’  net  Ilf  s.  it  is  \  necessary  to  use  a  very  finf.  qrid,  and  a 


corrf  spon.jnv;  incre&je  in  conputer  tme,  or  suffer  froM  reduced 
■I  c  r  I !  I  n  i.  V . 

'r,c  (1Pt^.ocls  described  here  enplov  an  irregular  triangular  grid,  instead 
of  a  regular  rectangular  grid,  so  it  is  possible  to  fit  ehpirical 
boundaries  with  high  precision  without  using  prohibitively  fine  gnus 
throughout  the  donain  of  the  solution.  The  penalty  for  this  is  a  slight 
increase  in  the  conple-.ity  of  the  algoi  ithn. 

Triangular  grids  have  been  extensively  used  for  elastic  and  plastic  flow 
problems  with  the  finite  elenent  Method  (see  for  exawple,  Strang  and 
Fix,  19;’3).  These  grids  have  also  been  used  in  Lagrangian  forMulations 
of  the  equations  of  hydrodynawics  by  Crowley  (1971),  Boris  et  al  (197S) 
and  Fritts  (1976),  but  the  Lagrangian  Method  iS  Most  suitable  for  flows 
where  the  total  deforMation  of  the  fluid  is  SMall.  Eulerian 
calculations  have  been  perforned  by  Sadoiirny  et  al  (1968),  UilliaMSon 
(  1968)  and  Thacker  ( 1 977) . 


This  report  describes  botti  ci  Me ttrodolo^v  tor  constructing  and  using 
ir  eTul:ir  tnanqiilar  irids  for  solving  ocean  flow  problens,  and  a  set 
0*  coNpLiter  progransfor  inplenenting,  testing  and  evaluating  the 
nethods.  The  tonplete  set  of  subroutines  is  given  later  as  an  appendix. 
Certain  details  of  tte  basic  algorithn  are  best  understood  tv  reference 
to  the  code  itself.  !he  code  has  been  designed  to  be  highly  nodcilar, 
so  that  sone  effort  «ust  be  expended  in  describing  the  interface  to  each 
nodule. 

To  n.i  inise  the  fU  ibilitv  of  the  code,  it  has  been  wrilten  in  a 
superset  of  fOMKftN  which  includes  nacre  expansions.  By  the  use  of  the 
file  inclusion  nacro,  tINLLUDE,  connon  block  naintenance  is  greatlv 
facilitated,  while  the  bulk  of  the  source  code  is  reduced.  Many 
conpilers  support  sone  forn  of  file  inclusion.  The  other  use  of  nacros 
in  the  code  is  to  support  paraneters.  Taraneters  are  used  for  values 
which  mist  take  the  forn  of  constants,  for  exanple  to  diMension  an 
array.  The  values  nav  need  to  be  changed  to  produce  different  versions 
of  the  code,  for  exanple  to  be  able  to  change  the  niinber  of  grid  points 
or  to  handle  different  hardware  configurations.  The  easiest  way  to  be 
able  to  enploy  nacros  is  to  use  a  sinple  pre-processor,  which  produces 
as  output  a  standard  FORTRAN  prograri,  which  anv  conpiler  can  then 


accept.  Alternatively, 


a  text  editor  nay  be  used  to  perforn  the  text 
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substitutions  inplie'j  bv  the  n^icros. 


Ill  the  work  thdl  i5  described  here,  a  sinpile  Macro  preprocessor  based  on 
that  of  Kerniqhan  and  Plauger  (1976)  was  used. 


The  following  predefined  Macros  are  assuMed: 


tHACROlNAHE, VALUE)  Define  new  Macro. 

t INCLUDE ( f  1  leiioMe )  Include  a  file  in  tlie 

source  code  at  this 
point. 

tDELTOK  Delete  the  next  input 

token,  norMally  a  new 
line  character. 


The  reMciining  Macros  are  defined 
uses: 


5I1AX 

Size 

IVHAX 

Size 

ICNAX 

Size 

SERAS 

Size 

in  the  code  and  have  the  following 

of  array  S. 
of  array  IV. 
of  array  IC. 
of  array  ERAS. 
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NAPMAXH 

IIAf'MAXV 

OUTLUN 

CHARACTER 


Hori:contJil  piige  size  in  chjiracter  widths. 
Vertical  page  size  in  character  hei'jhts. 
Logical  unit  nunber  for  the  output  device. 
Fype  for  character  datu«. 
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i  !HL  uKiD 

The  coMputat jonal  vjrid  coi.sif.ts  of  An  i rre^iil Ar  in An':)iilar  tesselation 
which  i;  tailored  to  approxi nate 1 v  correspond  to  the  shape  of  the  ocean 
hasiii  and  any  inlands  that  are  to  be  Modelled.  In  order  to  be  able  to 
CGiiveniently  handle  inch  a  grid  on  a  concuter,  we  mist  tornulate  a  data 
structure  that  represents  the  connect. vitv  and  the  -leoMetrv  of  the  gno. 
Fii-  =  t  let  us  in  trod  LI  re  two  definitions. 

A  yerte:;  is  defined  as  the  point  of  intersection  of  a  m.inber  of  grid 
lines. 

A  cell  IS  cl  triangular  region  bounded  by  three  grid  lines. 

Next  we  mist  ref  cr  specifically  to  the  code  (given  cSs  an  app>eTidix)  in 
order  to  See  r.ou  the  data  structure  i  org.inised.  Cell  connectivity 
1  nforiiation  is  field  in  the  an'av  1C  and  vertex  connectivity  infornation 
IS  held  in  the  array  IV. 

An  integer  expression  I  is  said  to  point  to  a  tell  if  1C(1)  is  the  first 
word  desrnt'ing  the  cell,  sinilarly  an  expression  J  is  said  to  pioint  to 
a  vertex  if  IV(J)  is  the  first  word  describing  the  vertex. 

In  addition,  a  nimber  of  words  of  real  storage  are  associated  uith  each 
cell  and  verle:  for  storing  physical  quaiitilies.  These  are  held  in  the 


it 
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arreiv  S. 


The  first  two  worijs  of  storage  for  vertices  sre  assimed  to  be  the 
coordinates  of  the  vertex.  Uhere  a  vector  field  is  to  be  represented, 
adjacent  words  in  S  are  used  to  hold  the  components  of  the  vector. 

As  an  example,  consider  the  shallow  water  equations.  The  physical 
storage  layout  associated  with  vertices  is  given  below: 

Uord 

0,1  Coordinates  of  the  verte'. 

2-8  Heights  for  lumping. 

9-15  Transport  weights  in  the  .x-di rection. 

16-22  Transport  weights  in  the  y-direction. 

23  Depth. 

24,25  Momentum. 

26-28  Momentum  flux. 

29-31  Lumped  depth  and  momentim  field  A. 

32-34  Lumped  depth  and  momentum  field  P. 

3j  32  Lumped  depth  and  momentum  field  C. 

38.39  Force. 

40,41  Coriolis  term. 

The  meanings  of  these  fields  are  described  elsewhere. 


k: 
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The  Icivout  of  connecti.ity  information  for  a  cell  is  as  -iiiven  below: 

1C(1)  Pointer  to  the  S  array. 

lhiI+1)  Fninter  to  the  tst  vertex  that  'iefines  the  cell. 

lC(I+2!  Pointer  to  the  2nd  ver((x  that  defines  th(  cell. 

IC(I+3)  Pointer  to  the  3rd  vertex  that  defines  the  cell. 


The  layout  of  connectivity  inforfiation  for  a  vertex  is  ^iven  below: 


U't  I ) 
IV(I  +  I ) 
IV(  1  +  2) 
1V( 1+3  ' 
IV(  1  +  4) 
IV(l+5) 
lV(I+3) 
IV( 1+7) 


Pointer  to  the  S  array. 

Nurtber  of  adjacent  vertices,  N. 

Pointer  to  the  1st  adjacent  vertex. 
Pointer  to  the  2iid  adjacent  vertex. 
Pointer  to  the  3rd  adjacent  vertex. 
Pointer  to  the  4th  adjacent  vertex. 
Pointer  to  the  5th  adjacent  vertex. 
Pointer  to  the  6th  adjacent  vertex,  or 
link  to  next  boundary  vertex. 


The  jrid  IS  assiined  to  have  the  sane  topoloqv  as  a  tessalation  of 
equilateral  triangles.  Arbitrary  boundaries  may  be  imposed  subject  to 
the  constraint  that  interior  boundaries  (islands)  do  not  have  acute 
anjles.  Since  the  jrid  is  deformed  before  use,  this  does  not  impose 
any  restriction  on  the  actual  shapes  of  islands  that  nay  be  handled, 
however  the  grid  would  be  more  than  usually  deformed  in  the  neighborhood 


J 
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of  such  fpcutp  f'ecilcirps. 

At  this  point  It  IS  iipijropr liite  to  exaniiie  portions  of  the  code  in 
detail.  ,0  we  will  examne  the  subroutines  that  plav  an  important  part 
in  constructinq  the  ^rid  data  structure. 

HEXTOF'  is  a  subroutine  that  constructs  a  grid  with  the  sane  topology  as 
a  hexagonal  region  divided  into  triangular  cells.  It  is  called  by: 

CALL  HEXTOP(N) 

Each  side  of  the  hexagon  is  divided  into  N  segnents. 

Bv  inspection,  we  «av  see  that  the  grid  has  6  N  triangular  cells, 
3(3N+1)N  edges  and  3(N+1)N+)  vertices. 

HEXrOP  generates  the  pointers  to  vertices  that  are  stored  in  IC,  by 
using  the  subroutines  GUC  and  GDC  which  generate  upward'  and  'downward' 
pointing  triangular  cells.  HEXIOP  counts  the  total  nunber  of  ceils 
allocated  in  NUMC,  and  the  nunber  of  vertices  in  NUMV. 

Next,  HEXTOP  generates  the  links  in  each  vertex  data  structure,  by 
exanining  each  cell  and  naking  sure  that  all  vertices  surrounding  a  cell 
are  joined  together  by  using  the  subroutine  VLINK. 


1 1 


Fin&llv,  HEArOP  dettTMViPS  which  vertu  .-5  lie  on  tiie  ijOLindury  of  the 
region  hv  co-.intici  the  ncMher  of  odjocent  vertices.  The  boiindfirv  points 
ore  coiitMined  in  o  linked  list  dots  structure. 

HEXTOP  IS  designed  ns  a  specific  tool  for  generating  a  class  of 
topologies  that  are  useful  for  investigating  the  properties  of  advection 
schenes  and  in  sinulations  of  anv  unbroken  region  that  has  approxinately 
circular  boundaries.  In  addition  HEXTOP  illustrates  the  nethods  that 
nay  be  used  to  generate  even  nore  conplicated  topologies,  such  as  grids 
with  inbedded  holes,  f ro«  the  hu«an  engineering  standpoint,  an 
interactive  conputer  svsten  with  a  display  screen  and  light  pen  would  be 
wore  convenient  for  generating  grids  fron  real  ocean  naps. 

GDC  and  GUC  are  cell  generating  subroutines  for  downward  and  up'ward 
pointing  cells  in  a  triangular  tesselation.  A  downward  pointing  cell  is 
generated  by  the  call; 


CALL  GDC(IP) 


where  IP  points  to  the  upper  left  hand  vertex.  It  is  assumed  that  the 
following  verte;..  IP+ItiSIZ.  is  the  upper  riqht  hand  verte;,  and  that 
the  Most  recently  defined  vertex  is  the  lower  verte:.  After  generating 
the  cell-vertex  links  for  a  new  cell,  GDC  returns  after  incre«enting  IP 


to  the  next  vertex. 


GUC  worlds  in  n  co«p leMentiiry  foilMoi!. 


ft  C  fi  i  1 

CALL  GUCdl'-J 

ftssurte?'  th?it  If'  points  to  the  upper  verte  ot  i<  cell.  Ihe  nost  recent 
verte  ;  is  sssuMed  co  te  the  lower  left  vertex,  so  GUC  generates  a  new 
vertei;  for  the  lower  right  hand  vertex  and  constructs  a  new  set  of 
cell -vertex  pointers.  IP  is  left  unchanged. 

Bv  alternating  calls  to  GUC  and  GDC  it  is  simple  to  construct  an 
arbitrary  tesselation  of  triangular  cells,  with  all  thie  reqtiired 
pointers  to  define  the  topology  of  the  grid. 

VLINK  15  a  subroutine  that  ensures  that  two  vertices  are  linked  together 
in  the  order  given.  It  is  called  as  follows; 

CALL  UL  INKUt'A.IVD) 

Uerte-.  It'A  is  examined  to  see  if  it  is  alrend'y  narked  as  having  lUB  as  a 
neighbor,  if  so  VLINK  exits  without  doing  anything  extra.  Otherwise 
the  vertex  count  for  lUA  is  increnenled  and  a  link  to  IVB  is  placed  in 
tf..  next  available  vertex  link  p  xsition.  ^o  correctly  link  two 


vertices,  two  calls  to  VLINK  are  required: 
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':''iLL  VI.  INKd'Jft.IVp  ! 
CALL  VLlNKlIVt.IV'.' 


The  subroutine  VSTO  iissigns  5lora'_:ie  spuce  to  the  v'-ftices  froi’i  the  array 
of  storage  S.  It  is  called  as  follows: 

CALL  VSTO(IORIG,NUIiS) 

The  first  word  allocated  is  ^(lORIG)  and  a  total  of  NlJnS  words  are 
assigned  for  each  verte:;.  NUnS  should  be  equal  to  the  ramber  of 
physical  variable  fields  required  in  the  solutions  of  the  equations  of 
hydrodynamics. 

CIRBND  IS  an  example  of  a  procedure  that  fixes  the  co-ordinates  of  the 
boundary  vertices,  in  this  case  by  distributing  them  uniformly  in  a 
circle.  It  may  be  called  without  any  arguments. 

CIRBND  IS  used  with  HEXTOP  and  SUGRID  t-.  define  a  circular  ocean  basin 
for  testing  purposes. 

An  assumption  is  made  in  CIRBND  that  all  boundary  vertices  lie  on  a 
single  exterior  boundary,  and  that  the  links  between  boundary  vertices 


follow  serially  around  the  boundary.  This  is  true  for  HEXTOP,  but  need 
not  be  so  for  other  topology  generation  subroutines.  A  more  general 
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boundrif-,'  fi:;in:t  routine  would  be  needed  tor  nultiple  boundaries,  and 
shouldprobotl V  be  enbedded  in  the  franewori  of  an  interactive  conputer 
avstert,  as  diacussed  elsewhere. 

VCPUN^  IS  a  subroutine  for  printing  out  the  links  associated  with  each 
verte::  and  cell.  It  is  called  without  anv  arguoents.  This  routine  has 
proved  useful  as  a  diagnostic  during  the  debugging  phase  of  new  topology 
generation  subrout i .its .  It  is  also  invoked  when  code  detects  an 
inconsistency  in  the  links,  caused  by  a  logic  error  in  the  topology 
routine,  or  none  usually  by  storage  corruption  caused  by  airciy  bound 
overflow  or  subroutine  argunent  inconsistencies,  which  traditioneil ly 
are  not  detected  by  FORTRAN  compilers. 

The  last  subroutine  in  the  grid  building  suite  is  SUGRID,  which 
calculates  the  co-ordinates  of  all  the  interior  vertices.  It  is  called 
as  follows: 


CALL  SUGRHKNITER) 

Interior  vertices  are  «oved  so  that  the  co-oiJi nates  of  each  vertex  is 
equal  to  the  average  of  those  of  all  its  neighbors.  This  involves 
solving  a  svsten  of  2N  linear  equations,  where  N  is  the  nunber  of 
interior  vertices. 


The  equations  are  anenable  to  a  relaxation  process  of  the  «ost 


strsi'jMforward  kind.  A  reliiXcilion  piirc'ineter  REITA  is  iised  to  inprove 
the  coiivei'sence  of  the  process. 


EftC'incol  tests  with  sinple  qnds  of  various  sices  have  shown  that  RELPA 
=  1.388  see«s  to  be  close  to  the  optinon  value. 

A  total  of  NITER  iterations  is  perforned,  and  for  each  iteration,  the 
naxinun  displacement  of  a  vertex  is  printed,  in  order  to  illustrate  the 
convergence  of  the  iterative  solution. 

On  the  basis  of  experience  with  other  relaxation  methO'js  applied  to 
elliptical  systems  of  partial  differential  equations,  it  is  believed 
that  SUGRID  is  unconditionally  stable  for  a  finite  range  of  RELPA, 

The  suite  of  subroutines  that  has  been  described  above  has  been 
specifically  written  lo  tonstruct  the  data  structure  for  a  circular 
ocean  basin,  but  simple  modification  to  CIRBNP,  for  example  would 
permit  irregular  quasi -c i > c u! ar  basins  to  be  treated.  If  the  basin  was 
grossly  dissimilar  to  a  circle,  or  had  a  different  topology,  for 
example,  containing  islands,  then  the  routine  HEXTOP  would  have  to  be 
replaced.  The  following  algorithm  is  proposed  for  the  [leneration  of 
grids  for  real -world  oceanographic  simulations: 

First  assemble  a  number  of  equilateral  triangles  to  approximatelv 
represent  the  region  of  interest.  For  example,  HEXIOP  uses  six 


t 
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triangles  in  the  ford  of  o  hexci^on  to  rep>resent  a  circle. 

The  ne>:t  step  li  to  subdivide  the  tritin-jles  into  a  tesselation  of 
srtoller  equilateral  triangles  until  it  is  estiwated  that  there  are 
sufficient  cells  to  resolve  the  solution  structure  that  is  desired. 

Triangular  cells  are  next  renoved  fron  the  interior  to  rei'itesent  islands 
and  froM  the  exterior  to  better  represent  the  ocean  shore.  The  vertices 
on  the  exterior  and  interior  boundaries  are  then  assigned  the  actual 
coordinates  of  associated  points  on  the  actual  coastline. 

The  final  step  is  to  relax  the  interior  points  to  obtain  a  grid  with  a 
srtooth  transition  of  cell  sire.  The  whole  procedure  would  best  be  done 
on  a  conputer  systen  under  interactive  connand.  A  video  graphics 
display  and  a  lightpen  would  be  dost  appropriate  for  adjusting  the 


coastline  points  in  order  nate  the  grid  as  uniforn  as  possible. 


4  THE  ALuOMIHlI  I  Ur<  HrcftOHrAMICS 

In  Host  fluid  dynamcs  siNulntion"  the  iiosl  troublesofle  terms  in  the 
equations  are  the  nonlinear  ones.  /he  adi/ection  of  a  scalar  field  may 
he  used  illustrate  the  alqorithrt  that  is  used  for  More  complicated 
cases. 


The  qoverninq  equation  for  a  scalar  field  is: 


c  V  .F  -  0 


Alternativelv,  hv  Stoles  theorem,  we  can  say  that  for  any  reqion  i, 
hounded  hv  a  curve  S,  we  have 

-  /' 

I  1 

♦  \  i-1^-  = 

n  ^  .  I- 

where  ^  is  normal  to  the  cm  ve  S,  and  is  the  average  value  of  ^ 

111  the  region  i . 

In  a  finite  representation  of  using  a  grid  of  points  on  a  triangular 
mesh,  it  IS  convenient  to  store  the  value  of  y)  at  the  vertices  of  the 
grid,  and  to  consider  a  region  surroundinq  each  vertex  which  we  will 
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col!  0  fhix  cell. 

A  flii,,  cell  15  an  irieqcilor  and  50Meti«e5  non-conve>:  polygon  that 
cont-ms  only  one  verte;:.  Each  side  of  the  polygon  starts  at  the 
^entoid  of  a  cell  adiacent  to  the  contained  vertex  and  finishes  on  the 
Mid-point  of  a  side  of  the  cell  that  passes  through  the  contained 
vertex.  It  is  clear  that  no  part  of  the  region  defined  by  the  grid  lies 
outside  of  all  flU/!  cells. 

A  flux  cell  surrounding  a  vertex  mth  N  neighboring  vertices  is  a 
polygon  with  2N  sides. 

Averages  of  a  variable  within  a  flux  cell  are  termed  '  UiMped'  values  and 
are  associated  with  the  contained  vertex.  If  the  values  of  the  variable 
are  to  be  represented  by  its  values  at  the  vertices  only,  as  is  usual 
for  anv  finite  difference  scheee,  then  we  nav  uniquely  approxinate  the 
variable  at  any  point  by  linear  interpolation  between  the  vertices  of 
the  cell  containing  the  point.  Uith  this  representation  of  the 
vaiiable,  we  nay  exactly  calculate  Iii«ped  values  bv  integrating  over 
the  flux  cell.  SiMilarly.  tfie  line  integral  above,  nay  be  exactly 
evaluated  around  the  flux  cell. 

Our  advection  equation  is  conservative  in  the  lutiped  values  pf  the 
scalar  variable.  Fluxes  are  derived  bv  first  evaluating  the  unlunped 
values  at  vertices,  evaluating  the  non-linear  flux  terns  at  each  vertex 


in'!  then  rert'orning  linear  interpolation. 


Thus,  we  have  a  sinple,  non-ant- i^uous  algorithm  for  discretising  the 
continuous  equations  of  hydrodvnanics. 

The  line  and  surface  integrals  of  interpolated  fields  nay  te  eiiactly 
representeo  by  a  weighted  sun  of  verfex  values  in  the  vicinity  of  the 
region  of  integration.  The  weights  do  not  change  with  tine,  so  the 
coNplicated  integrals  nay  te  siNplv  evaluated  with  a  NiniNun  of 
coNputati onal  effort. 

The  subroutine  CVLU  calculates  the  weights  A,  Bj  required  to  calculate 
the  vertex  lunped  values  Fi 

?i  ^  A  Fi  t  ^  F.i  Bi 

where  .i  is  a  vertex  which  neighbors  vertex  i. 

Consider  a  vertex  i  and  a  cell  k.  The  field  F  i.  define-j  at  I  and  the 
other  two  vertices  c  und  -i  such  that  the  cell  ^  is  defined  by  the 
triangle  if-q.  If  we  apcroxiNate  F  within  T  bv  linear  interpolation,  we 
Nay  integrate  F  withir  the  quadrilateral  ip'gq'  where  q  is  the  Nidpoint 
of  IP  and  p  IS  the  Nidpoint  of  ip.  and  g  is  the  centroid  of  ipq. 

This  integral  is 


■iitdaaBi 


vt*  =  .k  t  n/54  1-1  f  7/108  Fp  ivVIOB  Fgj 


where  K  is  the  area  of  cell  k. 

The  verteA  limped  value  of  F  at  i  is  the  sum  of  over  all 
quadrilaterals  as  above  which  have  i  as  a  vertex.  The  weiqhts  A  and  Bj 
«ay  thus  be  seen  to  be 


A  I  -11/54  y  .\k 

V  A 

B  |  =  /VI  08  ■  A. Ik 

t- 

where-^jk  =  A  k  if  j  is  a  vertex  or  k,  else  is.jk  -  0. 

CVLU  is  invoked  as  follows: 

CALL  C'JLUaPOS) 

IPOS  IS  a  pointer  to  a  set  of  7  variable  locations  which  are  to  contain 
the  values  of  A  and  Pi. 

The  subroutine  CTU  conputes  the  weiqhts^  and  Pj  to  calculate  the  line 
integral  T  of  a  vector  flux  F  around  a  flux  cell  surrounding  a  vertex  i 
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where 


and  j  15  a  vertex  neighboring  i. 

A  and  Bj  are  vectors  associated  with  vertex  i  and  are  stored  with  their 
x-conponents  at  IPOS  and  IPOS+J  respectively.  The  y-conponents  are 
stored  7  locations  after  the  corresponding  x-conponent. 

The  flux  cell  around  a  boundary  vertex  is  truncated  by  the  actual 
boundary,  so  that  part  of  the  line  integral  Must  be  evaluated  along 
sections  of  the  boundary,  however,  straightforward  use  of  linear 
interpolation  suffices  to  calculate  A  and  Bj.  The  details  of  the 
algorithn  need  not  be  described  here,  since  the  code  itself  serves  to 
define  the  Method. 

The  subroutine  is  called  as  shown; 

CALL  CTU(IPOS) 

where,  as  before,  IPOS  is  a  pointer  to  the  region  of  storage  that  is 
to  contain  the  weights.  In  this  case  H  words  are  required  to  store  the 


\ 


vector  weights. 
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The  correctness  of  the  cilriorithM  nay  be  tested  by  the  subroutine  TSTCTU, 
which  evaluates  the  line  integrals  for  a  flux  of  (1,0)  and  (0,1)  for  the 
boundcii'ies  of  the  flu;;  cells  contaiiunQ  each  vertex.  Uithin  the  Units 
of  rounding  errors,  the  two  integrals  are  zero. 

TSTCTU  IS  called  as  shown: 

CALL  TSTCTU(IPOS) 

where  IPOS  is  the  pointer  to  the  weights  generated  by  CTU. 

LUMP  IS  the  subroutine  that  calculates  the  integral  of  a  physical 
variable  field  over  the  flux  cell  that  surrounds  each  vertex. 

The  subroutine  LUMP  is  invoked  as  follows; 

CALL  LUMP( IP0S,IFR0M,IT0) 

IPOS  is  a  pointer  to  the  set  of  weights  that  have  previously  been 
calculated  by  CVLW,  and  IFROM  and  ITO  are  pointers  to  the  source  and 
destination  fields  respectively. 

The  lunped  variable  is  stored  at  location  ITO  after  calculating 
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where  A  and  B i  are  the  weights  and  I'j  is  the  value  of  the  field  at 
verte;i  .i.  j  is  a  neighbor  verte>;  of  vertex  i. 

UNLUMF'  IS  a  subroutine  which  is  the  fornal  inverse  of  LUHF'.  It  will 
recover  the  original  values  of  a  field  that  has  previously  been  smoothed 
by  LUMP.  It  15  invoked  as  shown: 

CALL  UNLUMP(IPOS,IFROM,ITO,ERR) 

IPOS,  IFROM  and  ITO  have  the  sane  weaning  as  in  the  above  description 
of  the  subroutine  LUMP. 

The  systew  of  linear  equations  for  Fi,  given  Fi  is  solved  by  relaxation 
Methods,  with  a  final  accuracy  estinated  to  be  on  the  order  of  ERR  if 
possible,  otherwise  an  error  Message  is  printed. 


Ihe  subroutine  ADUECT  is  called  as  follows: 

CALL  ADVECK IPCN.IPCL,IPC,IPU,IPU,DT  ) 

This  subroutine  perforws  the  basic  advection  of  a  limped  scalar  field  at 
location  IPCL  using  the  unliiwped  velocity  whose  x-  and  y-co«pcnents  are 
located  at  IPU  and  IPUH.  Transport  weights  are  located  at  IPU.  The 


tiwe  interval  for  integration  is  BT  and  the  new  scalar  field  is  returned 


to  the  location  IPCN. 


The  subroutine  ADVECT  modifies  a  lufiped  liOMentun  field  to  ensure  that 
the  corresponding  nonentuN  values  on  the  grid  boundaries  are  zero.  This 
15  equivalent  to  inposing  a  rigid  boundary  condition  at  the  edge  of  the 
conputational  donain. 

The  subroutine  STEP  is  called  as  follous: 

CALL  STEP(IPLN,IFULN,II-'VLN,IPLO,IPULO,IPVLO,IPL, 

IPUL,IPVL,IP,IPU,1PV,IPUU,IPUV,IPVV, 

IPFX,IPFY,IPFXL,IPFYL,IPU,DT,F,CF, 

USX.USY) 

The  STEP  subroutine  advances  the  fields  of  the  physical  variables  one 
tine  step  for  the  shallow  water  equations.  The  values  of  the  dependent 
variables,  height  and  horizontal  nonentun  conponents,  at  the  start  of 
the  step  are  denoted  with  a  suffix  '0'  for  Old  while  the  values 
corresponding  to  the  end  of  the  step  have  a  suffix  'N'  for  New".  STEP 
first  unlunps  the  fields  in  order  to  calculate  pressure  and  Reynolds 
stress  terns,  then  uses  the  transport  algorithn  and  finally  inposes  the 
boundary  conditions. 


The  subroutine  TRNSPT,  which  is  called  as  shown; 
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CALL  TRNSPT(INEU,IOLIi,lFX,IFY,  IPU.IiT  ) 

evrtlustes  the  changes  to  a  lunped  field  for  a  tine  step  DT,  The 
original  luMped  field  is  located  at  lOLP.  The  field  after  Modification 
for  the  effects  of  transportation  is  returned  to  location  INEU.  The 
transportation  is  effected  by  a  flux  field  whose  x-  and  y-conponents  are 
located  at  IFX  and  IFY  respectively. 

The  shallow  water  equations  contain  several  terMs  which  are  not  of  an 
advective  nature.  These  terns  are  generally  easier  to  treat  in  a 
nunerical  schene  than  the  non-linear  advective  terns.  Ue  will  refer  to 
the  non-advective  terns  as  forcing  terns.  The  subroutine  FRF  includes 
the  effects  of  the  following  forcing  terns: 

1)  Frictional  tern  with  a  coefficient  of  friction 
Cf  and  a  linear  dependency  with  the  fluid 
speed. 

2)  Coriolis  tern  with  a  paraneter  F. 

3)  An  externally  inposed  stress  (USX.USY)  which 
nay  be  regarded  as  representing  the  effects  of 
wind  on  the  fluid. 


I 

1 


MiMiiiiiiili 


The  subroutine  c<'illed  witli  the  rirqudents 


FRFdF'L.IPUL.IPVl  .IP,IPU,IPV.ir-'FX,IPFY,IF'FXL, 
lPFyL.DT,F,CF,USX,USY) 

IP,  IPU  iind  IPy  nre  pointers  to  unliiMped  depth  and  nofientuM  components. 
IPl,  IPUL  and  IPVL  are  the  locations  of  the  corre5pondin3  lunped 
fields.  IPFX  and  IPFY  ai e  fields  used  to  store  the  conbined  frictional 
and  external  stresses.  IPFXL  and  IPFYL  are  used  to  locate  tfie  limped 
values  of  the  previous  two  fields.  Upon  exit  fron  FRF  the  momentim 
fields  are  Modified  to  reflect  the  ti«e  integrated  effects  of  the 
forcing  terns  over  a  tine  interval  DT. 

Two  Main  progran  are  included  in  the  appendix.  The  first  is  used  in 
simple  test  solutions  for  the  color  equation,  in  order  to  excmine  the 
stability  and  accuracy  of  the  basic  advection  algorithm. 

The  second  subroutine  is  a  driver  for  the  shallow  water  equations  and 
includes  calls  to  the  graphics  routines  which  are  described  in  the  next 
section.  The  subprogran  Makes  two  calls  to  STEP,  since  the  basic  tine 
differencing  schene  that  is  enployed  is  the  robust  pseudo-backward  Euler 


Method . 
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'i  our  PUT  ROOIINES 

In  debU'Vjin-j  the  rode  and  evsluatin-q  the  perf ornanc e  of  the  niifierical 
schewe  that  it  embodies,  ^raohical  output  is  an  invaluable  tool. 
SiNple  listings  of  field  values  give  little  indication  of  the  actual 
behavior  of  the  tiodel  being  sinulated,  especially  when  an  irregular 
grid  15  used.  A  nuMber  of  routines  have  been  produced,  which  serve  in 
thenselves  as  valuable  aids,  while  at  the  sane  tine  forning  a  set  of 
prinitives  for  none  elaborate  displays. 

The  QPRCQN  subroutine  is  a  ample  display  procedure  for  producing 
contour  plots  of  a  field  on  a  hardware  device  such  as  a  lineprinter.  It 
IS  called  as  follows: 


CALL  QPRCON(IFLB,XMlN.XhAX) 


IFLD  15  a  pointer  to  the  field  that  is  to  be  dispUayed.  A  single 
character,  either  a  digit  or  a  plus  or  ninus  sign  is  placed  on  the  page 
at  a  position  corresponding  to  a  vertex  location.  The  character 
represents  the  value  of  the  specified  field,  with  0  corresponding  to 
XMIN  and  9  corresponding  to  XhAX.  Internediate  values  are  represented 
by  an  appropriate  decinal  value  using  a  linear  transfer  function. 
Values  less  that  XttIN  are  shown  by  a  Minus  sign  while  values  larger  that 


XHAX 

are  shown 

with 

a  plus  sign. 

If  the  grid  is 

coarse,  the  display 

will 

be  sparse. 

but 

QPRCON  has  the 

advantage  of 

representing  both  the 

grid  structure  and  the  field,  so  that  ludgenents  nay  be  «ade  as  to  the 
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reality  of  features  with  scales  near  that  of  the  basic  qri'j. 

PRCON  IS  a  subroutine  that  is  functionally  sihilar  to  OPRCON,  having 
ihe  sane  argunent  list,  e/cept  that  the  field  values  for  space  between 
vertices  are  also  displayed.  Each  character  space  on  the  output  page  is 
napped  into  the  grid,  and  if  the  point  lies  within  the  boundary  of  the 
grid,  the  value  of  the  input  field,  IFLD,  is  evaluated  by  linear 
interpolation.  The  running  tine  for  PRCON  nay  be  nurh  larger  than  that 
for  QPRCON,  especially  when  the  the  nunber  of  vertices  is  less  that  the 
nunber  of  resolvable  positions  on  the  output  page. 

The  routine  CONPLT  generates  a  contour  line  plot  of  the  specified  field. 
It  IS  called  by  neans  of: 

CALL  CONPLT(IFLD,CONVAL,N) 

Contour  values  are  printed  for  values  CONVALd),  1  =  1  ,N.  CONPLT  is 
suitable  for  all  types  of  contour  plotting  devices  such  as  pen  plotters 
and  electrostatic  printer-plotters,  since  it  nakes  use  of  only  one 
svsten  dependent  subroutine,  LINE,  for  drawing  a  straight  line 
segnent . 

The  DRAUV  subroutine  is  called  as  follows: 


DRAUV(IPL,IPUL,rPVL,IP,IPU,IFV,ASCALE,IPU) 
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The  IifiAUV  subroutine  9enfer-jtes  ;i  display  of  the  velocity  field  iwplied 
by  the  lunped  depth  and  Mo«entu«  fields  located  at  IPL,  IPUL  and  IPVL. 
The  fields  IP,  IPU  and  IPV  are  used  as  workspace  for  calciilatin-j 
unlunped  values  of  the  physical  variables. 

The  COMMON  variable  SCALE  is  used  to  transforn  froe  ^rid  coordinates  to 
plotter  coordinates.  After  unlunping  the  fields,  the  velocity 
cofiponents  are  ei^tracted  fron  the  Monentun  fields.  At  each  vertex  of 
the  grid,  an  arrow  is  drawn  with  a  length  proportional  (using  the  scale 
factor  ASCALE)  to  the  velocity  at  the  specified  point.  No  attenpt  is 
Made  to  draw  arrows  that  are  so  snail  as  to  be  unresol vable. 

DRGRID  IS  a  subroutine  that  has  no  argiments  but  that  nay  be  used  for 
drawing  the  triangular  grid  on  a  scale  sinilar  to  the  other  plotting 
displays.  so  that  it  nay  be  used  to  overlay  contour  or  vector  field 
representations. 

Fho  subroutine  [iRCELL  is  sinilar  to  DRGRID,  also  having  no  argiments, 
e.cept  that  it  draws  the  outline  of  the  flux  cells  that  are  used  for  the 
line  integrals  in  Stokes'  theoren. 

EUPLOr  IS  a  subroutine  with  the  following  calling  segnence: 


CALL  EUPLOT(IP.ZMAX,ZNIN) 
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The  EUPLOT  sutroiitine  is  used  to  generiite  a  plot  of  the  field  located  at 
IP,  along  the  line  T  =  0.  The  actual  plot  15  produced  by  printing 
characters,  in  a  way  suitable  for  use  with  a  linecrinter.  The  rta;;iMUM 
and  niniMUM  values  that  nay  be  displayed  on  the  page  are  given  by  ZMAX 
and  ZMIN  respectively. 

INTPOL  IS  the  basic  interpolation  subroutine  used  by  the  other  graphics 
subroutines  for  extracting  field  values  fron  arbiti'ary  positions  on  the 
grid  bv  neans  of  linear  interpolation.  It  is  called  as  follows: 

CALL  INTPOL(IFLli,X,Y,Z,FAIL) 

The  postion  specified  is  given  by  the  coordinates  (X,r)  and  the  value  of 
the  field  indicated  by  IFL  is  returned  in  the  variable  Z.  If  the  point 
lies  outside  the  boundarv  of  the  grid,  the  logical  variable  FAIL  is  set 
to  betrue.  INTPOL  works  by  exanining  each  cell  in  turn  until  one  is 
found  that  contains  the  required  point,  so  it  is  relatively  tine 
consuning,  especially  when  it  is  called  nanv  tines,  as  in  the 
subroutine  PRCON. 

The  routine  could  be  increased  in  speed  by  introdcicing  a  relatively 
coarse  rectangular  grid  and  associating  a  linked  list  of  triangular  cell 
inmbers  with  each  rectangular  cell.  The  search  could  then  be  perforned 


on  a  nuch  shorter  list,  the  particular  list  being  determined  by  direct 


r 
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conputation  frort  X  and  r.  This  would  be  anal090us  to  hash  coded 
searches  in  standard  table  look  up  problens,  with  the  added  advantage 
that  each  list  would  be  approxinately  the  sane  size,  so  the  look  up 
tine  would  be  both  short  and  predictable. 

INTRI  IS  function  that  return  a  value  of  LOGICAL  type.  It  «ay  be 
referenced  by; 


IMTR1(X,Y,ICELL) 

The  returned  value  is  only  true  if  the  point  (X,Y)  lies  inside,  or  on 
the  boundary  of  the  triangular  cell  nunber  ICELL. 

The  subroutine  hAXHIN,  which  is  called  by: 

CALL  MAXNIN(IFLD) 

prints  out  the  naxinun  and  nininun  values  of  the  field  specified  by 
IFLD, 


The  function  NEIGH  returns  a  LOGICAL  value  if  the  two  vertices  specified 
by  IP  and  10  are  neighbors  when  it  is  called  by  the  sequence: 


NEIGH(IP,IQj 
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The  subroutine  PRFLDS,  called  by 

CALL  PRFLDSdl  ,12) 

prints  out  the  values  of  all  the  fields  in  the  range  froe  II  to  12  in 
tabular  forn  for  each  vertex.  If  needed  the  table  is  split  into  several 
parts  so  as  not  to  exceed  the  uidth  of  the  lineprinter  page. 

The  VCDUMP  subroutine  displays.  on  the  printer,  all  the  links 
associated  with  the  cell  and  vertex  data  structures.  VCDUHP  does  not 
take  any  argunents. 

CONSRV  IS  a  subroutine  called  as  follows: 

CALL  CONSRVdPROff) 

This  diagnostic  routine  evaluates  and  prints  the  sun  of  the  values  of  a 
variable  at  each  vertex  of  the  grid.  If  the  variable  is  a  lunped  field, 
the  sun  is  the  spatial  integral  of  the  corresponding  unlunped  field. 
The  routine  has  proved  itself  useful  in  code  testing  by  evaluating  the 
integral  of  fields  that  should  be  fornally  conserved  in  tine  by  the 
systen  of  equations. 

The  COPY  utility,  which  is  called  by: 
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CALL  COPY(ITO,1FROI1) 

is  used  in  nany  places  throughout  the  code  to  transfer  a  physical 
variable  field  froH  one  location  IFROH  to  another  location  ITO. 

The  subroutine  FLDAV  is  a  utility  procedure  which  is  called  by: 

CALL  FLDAV(IPOSA,IPOSB,IPOSC) 

It  calculates  the  average  of  two  fields  indicated  by  IPOSB  and  IPOSC  and 
returns  the  result  to  the  location  IPOSA. 

INITF  is  a  subroutine  used  to  set  up  initial  fields  of  unit  depth  and 
zero  nonentufl.  It  is  called  by: 

CALL  INITF(IP,IPU,IPV) 

The  subroutine  SBROT  imposes  a  velocity  field  appropriate  for  a  solid 
body  rotation  with  an  angular  velocity  of  0HE6A.  The  resulting 
conponents  are  stored  at  locations  IPOSU  and  IPOSU+1  when  it  is  called 
by : 


CALL  SBROTtOHEGA, IPOSU) 


The  TFIELD  subroutine  sets  up  a  scalar  test  field  uhich  is  a  conpact, 


I 
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continuous  and  differentiable  piecewise  bivariate  cubic,  in  the  shape 
of  an  off-center  'hunp'.  This  field  has  proved  useful  for  studyin9  the 
advection  of  a  scalar  field  under  solid  body  rotation.  Although  this 
problen  is  analytically  trivial  since  the  field  nerely  rotates  with  the 
flow  without  changing  shape,  it  is  quite  non-trivial  conputationally 
and  provides  a  powerful  test  for  conparing  the  efficacy  of  different 
nunerical  advection  schenes.  It  nay  be  invoked  by  the  call: 


CALL  TFIELD(IPOS) 
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6  RESULTS 

A  sequence  of  nuMencal  experiMents  uere  perforned  with  the  code  durin'j 
its  develop«ent  in  order  to  provide  answers  to  the  followinq  questions: 

1)  Can  a  sinple  al^orithn  be  developed  which  is 
capable  of  solving  the  equations  of  hydrodynanics 
with  the  physical  variables  approxinated  by  values 
on  an  irregular  grid? 

2)  Uhat  are  the  Units  of  the  stability  of  the 
nunerical  schene? 

3)  Uhat  is  the  accuracy  of  the  schene,  especially  in 
the  treatnent  of  the  doninant  non-linear  advective 
terns? 

4)  Uhat  conputing  resources  are  required  to  use  the 
algorithn? 

5)  How  does  the  algorithn  conpare  to  orthodox  finite 
difference  schenes’ 

The  answer  to  question  I)  is  clearly  yes.  The  schene  that  has  been 
described  is  certainly  sinple  and  elegant.  Sone  slight  conpl ications 


are  inherent  in  using  an  irregular  grid,  but  all  the  geonetrical  and 


topological  factors  are  resolved  once  at  the  start  of  a  siMulation  and 
are  incorporated  in  the  weight  terns  and  the  cell  and  vertex  links. 
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During  execution  the  schene  is  fully  explicit,  except  for  the  sinple 
relaxation  schene  that  is  used  for  unlunping'  the  physical  variable 
fields . 

One  especially  satisfactory  feature  of  the  algorithn  is  that,  for  the 
treatnent  of  spatial  derivatives  at  least,  there  are  no  arbitrary 
decisions  to  be  nade  in  setting  up  the  difference  equations.  Host 
finite  difference  schenes  are  overdeternined,  in  the  sense  that  several 
different  representations  have  sinilar  spatial  order  of  accuracy,  with 
often  no  obvious  neans  of  resolution.  Also  the  present  algorithn  nay  be 
readily  generalised  to  nore  conplicated  situations. 

The  answer  to  question  2)  is  not  capable  of  analytical  solution,  so  we 
nust  perforn  experinents  and  nake  cautious  inferences.  Host  finite 
difference  schones  that  are  explicit  in  tine  differencing  nust  have  sone 
constraint  on  the  naxinun  tine  step  in  order  to  satisfy  stability 
criteria.  If  these  stability  criteria  are  not  satisfied  then  a 
conputational  node'  solution  nay  be  obtained.  These  frequently  exhibit 
high  spatial  frequency  oscillations  and  bear  no  resenblance  to  the 
solutions  of  the  underlying  differential  equations. 


Sone  nunerical  schenes,  such  as  the  Dufort-Frankel  schene,  while  being 
explicit  in  tine,  have  no  fornal  constraint  on  the  tine  step.  However 


they  achieve  stribilitv  at  the  expense  of  introdiic  in q  increasin'^  anotmts 
of  artificial  diffiisivity  as  the  tine  step  is  nade  larger.  In  addition, 
even  schenes  that  nake  use  of  inplicit  tine  differencin9,  such  as  ADI 
nelhods  are  still  in  practice  linited  by  the  need  to  naintain  adequate 
accuracy  in  the  tine  integration. 

A  nimber  of  nunerical  tests  have  been  perforned  both  for  the  scalar 
advection  code  and  the  shallow  water  equations.  For  seni -regular  grids 
the  naxinun  tine  step  for  stability  was  obse)'ved  to  be  inversely 
proportional  to  the  nean  inter-vertex  spacing,  while  for  nore  irregular 
grids  the  naxinun  stable  tine  step  was  found  to  be  none  closely  related 
to  the  inverse  of  the  nininun  inter-vertex  spacing.  The  latter  result 
is  hardly  suprising  but  it  does  nean  that  when  constructing  an  irregular 
grid  to  follow  a  coastline  of  creat  conplexity,  it  pays  not  to  use 
gratuitously  snail  cells.  However  if  the  boundary  points  are  reasonably 
distributed  along  the  coast,  and  the  grid  topology  is  appropriate  to 
the  conputational  c'onain  the  use  of  SUGRID  or  its  generalisation  will 
ensure  that  no  excessively  snail  inter-vertex  spacings  are  generated. 

Overall,  it  has  been  found  that  the  stability  behaviour  of  the 
triangular  grid  schene  is  qualitatively  and  quantitatively  sinilar  to 
nore  conventional  rectangular  grid  finite  difference  schenes. 

Sone  answers  to  question  3)  have  been  obtained  by  solving  the  color 
equation  for  a  solid  body  velocity  field.  Such  an  equation  has  a 
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trivuil  analvticsJ  soii.itioo  so  wr  nay  investigate  the  accuracy  of  the 
scheme  in  its  irealBent  of  the  all  inportant  advection  terN,  For  a 
circular  hasin  of  unit  radius,  a  solid  hody  rotation  law  and  an  initial 
field  deternined  hv  the  sutroutine  1FIELD,  'jood  results  were  obtained 
even  with  a  very  coarse  qrid.  For  exanple  with  a  typical  inter-vertex 
spacing  of  .125,  the  bell  shaped  test  field  only  suffered  about  102 
error  after  a  conplete  circuit  of  the  grid.  Since  positivity  is  not 
ensured  bv  the  schene,  so«e  snail  negative  and  positive  ripples  were 
observed,  especially  in  the  wake  of  the  “hunp",  but  these  always 
renamed  snail  in  nagnitude.  The  algorithn  has  second  order  spatial 
accuracy  so  in  the  linit  of  fine  grids,  error  estinates  nay  be  easily 
calculated. 

Question  4  nay  be  answered  both  on  theoretical  estinates  and  on  the 
results  of  the  nunerical  experinents  that  have  been  perforned.  If  we 
take  as  a  neasure  of  the  grid  conplexity  sone  nunber  N,  for  exinple 
NUhHIV  as  used  in  the  code,  or  the  reciprocal  of  the  inter-vertex 
spacing.  Then  the  nunber  of  grid  points  for  which  calculations  must  be 
perforned  is  N»»2.  The  solution  of  the  systen  of  linear  equations  for 
the  unlunping  operation  requires  on  the  order  of  Nf*3  operations,  since 
It  IS  iterative,  and  infornation  nust  be  propagated  over  the  whole  grid 
in  order  to  effect  a  solution.  For  each  tine  step,  the  conputational 
effort  has  two  factors,  one  proportional  to  the  square  of  N  and  the 
other  proportional  to  the  cube  of  N.  For  noderate  values  of  N  up  to 


about  20, 


the  first  tern  is  doninant 


For  fine  scale  sinulations, 


however,  the  effort  of  the  UNLUMP  operation  will  iJo«inate  the  solution. 
The  naxiNun  pernitted  tine  step  vanes  inversely  with  N,  so  a  conplete 
sinulation  will  take  a  tine  that  is  proportional  to  the  cube  of  N  for 
coarse  'jiricts  and  the  fourth  power  for  sufficiently  fine  grids. 

There  is  no  clear  cut  answer  to  the  final  question.  Sone  cartesian  grid 
Methods  Maintain  H**2  tining,  and  so  for  very  fine  grids,  will  be  More 
efficient  than  the  irregular  triangular  grid  Method.  On  the  other  hand, 
irregular  boundaries  May  require  excessively  fine  grids  for  the 
cartesian  grid  Methods,  so  their  theoretical  edge  May  not  be  Maintained 
in  practice.  It  has  been  denonstrated  that  irregular  triangular  grid 
Methods  can  be  applied  to  the  equations  of  ocean  flow,  without  undue 
penalties  in  coMputer  resources,  or  conplexity  of  the  final  code. 
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Append  i>: 


Typical  set  of  nncro  definitions 

tHACR0;',t$DELT0K3)" 

$HACR0(SHAX,6510)‘ 

tnACR0<IVHAX,173A)' 

$«ACR0(ICrtAX,1668)' 

»HACR0(HAPHAXH,72)“ 

♦  NACR0(hAP«AXV,22)'^ 
(HACRO(SERAS,n‘ 
tHACR0(0UTLUN,5)' 
tHACRO(INLUN,S)~ 
t«ACRO(CHARACfEfi. LOGICAL  ♦!)' 
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flam  progran  for  driving  the  scalar  advection  code 


$INCLUDE(CH(lC.MAC)»Ii£LTOK 
TYPE  1 

1  F0R«AT(  '  N,DT-') 

ACCEPT  »,M,0T 
URlTE(aUTLUN,2)  N,0T 

2  FORNATf'  N,DT',I4,F12.5) 

CALL  CTEST(N,DT) 

STOP 

END 

SUBROUTINE  CTEST(NUhDIV,DT) 

URITE<1,96)  NUNDIV.DT 
96  FORMATC'  N ,  DT  ' ,  I6,F M . 6 ) 

0HEGA=6. 263185 
NSTEPS=1 ./DT 
CALL  HEXTOP(NUNDIV) 

CALL  VST0(I,30) 

CALL  CIRBND 
CALL  SUGRID(20) 

IPLW=2 

CALL  CVLU(IPLU) 

IPU=9 

CALL  CTU(IPU) 

C  CALL  TSTCTU(IPU) 

CALL  GRDLEN 

IP  =  23 

IUs24 

IV=25 

IPL=26 

IPN»27 

IUL=28 

IVL=29 

CALL  TFIELD(IP) 

CALL  SBROTfOHEGA.IU) 

CALL  LUNPdPLM.lF  .IPL) 

C  Remove  the  Cs  free  the  next  three  lines  to  inpose  rigid  b.c. 

C  CALL  LUHP<IPLU,IU,IUL) 

C  CALL  LU«P(IPLU,IV,IVL) 

C  CALL  RIGBND<IUL,IVL,IU,IV) 

DO  to  I=I,NSTEPS 

CALL  AOVECT(IPN,IPL,IP,IU,IPU,BT) 

CALL  COPy(IPL,IPN) 

CALL  UNLUHP<IPLU,IPL,IP,1 .E-5) 

IF  (NOOdJO)  .EO.  0)  CALL  HAXHINdP) 

10  CONTINUE 

CALL  NAXNINdP) 

RETURN 

END 
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SUBROUTINE  6RDLEN 
$INCLUDE(TRG.DCL)' 

J=1 

NSEC=0 

XHIN>1.EtO 

XNAX:=0. 

XfiEAN=0. 

DO  1  i-i,Nunv 

JP=J+1 

NV=IV<JP) 

JA=IV(J) 

X=S(JA) 

Y=S(JA+1 ) 

DO  2  K=I,NV 

JK=JP+K 

JC=IV(JK) 

JB=IV( JC) 

XX=S<JB) 

YY=S(JB+I ) 

D=SQRT( (X-XX)**2t(Y-YY)**2) 

NSEG^NSEG^I 

XHEAN=XHEAN^D 

IF  <0  .GT.  XNAX)  XnAX>D 

IF  (D  .LT.  XHIN)  XHIN>0 

2  CONTINUE 

1  J=J+IVSI2 
XNEAN»XNEAN/NSEG 

3  FORNAK'  HIN,HEAN,MAX',3Ft2.5) 
URITE(1,3)  XNIN,XNEAN,XHAX 
RETURN 

END 

SUBROUTINE  TFIELDdPOS) 
«INCLUDE(TRG.DCL)' 

YC*0. 

XC  =  1. 

RC=.9 
1  =  1 

DO  1  ii>i,Nunv 
K=IV<I) 

X=S(K) 

Y=S(K  +  1 ) 

R=S8RT((X-XC)*»2+(Y-YC)**2)/RC 

C=0. 

IF  <R  .LT.  1.)  C=(2.*R-3.)*R*R*1. 
S(K+IPOS)*C 
1  I=I+IVSIZ 
RETURN 
END 


1 


I 


Nam  proqrcin  far  the  shailov  waiter  equations  code 


tlNCLUDEtTnAC.HAOtHELTQK 
DIMENSION  ERAS(SERAS) 
DIMENSION  CONOAL(IO) 
LOGICAL  PLOT 
PL0T=. FALSE. 

ASCALE=1 . 

DT=.04 

CALL  INPUK'DT _ ',DT) 

X  =  2 

CALL  INPUK'NUMDIV  ,X) 

NUN0IV=X 

X=101 

CALL  INPUT!  NSTEPS'jX) 
NSTEPS=X 
X  =  10 

CALL  INPUTCNCONT  '  ,X) 
NC0NT=X 
DO  2  1=1,10 
2  C0NVAL(I)=I«.2-.1 
CALL  HEXTOP(NUHDIV) 

CALL  VST0<1,t2) 

CALL  CIRBND 
CALL  SUGRID<20) 

IPLU*2 

CALL  CVLU(IPLU) 

IPU=? 

CALL  CTU(IPU) 

CALL  TSTCTU(IPU) 

IP  =  23 

IPU=24 

IPV=25 

IPUU=26 

IPUV=27 

IPVM=20 

IPLA=2? 

IPULA=30 

IPVLA=31 

IPLB*32 

IPULB=33 

IPVLB=34 

IPLC=35 

IPULC=3A 

IPWLC=37 

IPFX=38 

IPFY=3? 

IPFXL=AO 

IPFYL=41 

CALL  IMITF<IP,IPU,IPV) 
CALL  LUHP(2,IP,IPLA) 
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CrtLl  LUrtPCJ.IF'UjlF'ULA) 

CALL  LUMP(2,IPV,IPVLA) 

CF=0. 

CF  =  1  . 

CALL  INPUff  CF - .CF) 

F^O. 

CALL  input;  F .  ,F) 

USX^.3 

usy=o. 

CALL  input;  USX...  .USX) 

CALL  input;  UST...  ,USV) 

IF  ;PLOT)  CALL  PLO TS ; ERAS , 900 ,5. 5 ) 

DO  10  I=1,NSTEFS 

IF  («OD(I-l  .NCOfll  )  .NE.  0)  GO  TO  9 
CALL  UNLUnP;2.IPLA.IP.I .E-5) 

CALL  QPRCONdP.O.,:.) 

CALL  EUPL0T(IF,2..0.) 

IF  ;PLOT>  CALL  0RIGIN(5,5,0. ) 

IF  <PLOI)  CALL  drauv;ipla,ipula,ifvla,ip,ipu,ipv.ascale.iplu) 
9  CONTINUE 

CALL  step;  1FLB,IPU1.B,IPVLB, 

1  IPLA,IPULA,IPVLA, 

2  1PLA,IPULA,IPVLA, 

3  IP,IPU,IPV.IPUU,IPUV,IPVU, 

4  IPFX,IPFr,IPFXL,IPFYL,IPU,DT,F,CF,USX,USY) 

CALL  NAXNINdP) 

CALL  NAXMINdPU) 

CALL  NAXMINdPV) 

CALL  STEPdPLC,IPULC,IPVLC, 

1  IPLA,IPULA,IPVLA, 

2  IPLB,IPULP,IPVLB. 

3  IP,IPU,IPV,IPUU,IPUV,IPVV, 

4  IPFX,IFFY,IPFXL. IPFrL.IPU,Dr,F,CF,USX,USY> 

CALL  COPYdPLA,IPLC) 

CALL  COPYdPULA.IPULC) 

CALL  COPYdPVLA.IPVLC) 

10  CONTINUE 

IF  (PLOT)  CALL  ORIGINS. 5. 0. ) 

IF  (PLOT)  CALL  ENDPLT 

STOP 

END 

SUBROUTINE  INPUT ( STRING ,X  ) 

CHARACTER  3TRING<6) 

URITE(0UTLUN,1 )  STRING, X 

1  format;  Enter  value  for  ',AA1,  Perhaps' , 1 F 1 2. 3 ) 

ACCEPT  *,X 
RETURN 
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Connon  block  TRG.DCL 


COHHON  /TRG/  S< SHAX) , IV( IVHAX ) . IG ( ICHAX ) , 
.  SCALE, 

.  NEUV,NEUC,IVSlZ,ICSIZ,NUHC,NUnb,NBV,IBVF 


Library  of  subroutines  for  triangular  9rid  codes 


$INCLUDE(THAC.HAC)tDELrOK 

SUBROUTINE  ADVECT(IPCN,IPCL,IPC,IPU,IPU,DT) 
$IMCLUDE(TRG.DCL)‘ 

J=1 

DO  1  I=1,NU«U 
K=IV( J) 

MV=IV(J+1 ) 

KIPC=K+IPC 

KIPU=K+IPU 

KIPU=K+lPg 

KIPCN=K+IPCN 

S(KIPCN)=S(KIPC)»(S(KIPU)  t'S(KIPU)*S(KIPUH)»S(KlFU+7>) 
DO  2  IN=1,NV 

JP=IV(JIN  +  t ) 

KP=IV( JP) 

K1PCN=K+IPCN 

KPIPC=KP4IPC 

KPIPU*KP+IPU 

KIPUIN-K4IPU4IN 

2  S<KIPCN)sS(KIPCN)4S(KPIPC)*(S<KPIPU>*S(KIPM1N)* 

1  S(KPIPU4t )«S<KIPUIN47n 
KIPCN=K+IPCM 
KIPCL=KtIPCL 

S<KIPC«)»S<KIPCL)-DT*S(KIPC«) 

I  J=J4IVSIZ 
RETUR- 
END 


SUBROUTINE  CIRBND 
$INCLUDE(TRG.DCL)“ 

DTU=6.283184/NBV 

J*IBVF 

TH=0, 

DO  t  1=1, NBV 

IF  (J  .EO.  0)  GO  TO  ? 

IF  (IV(J41)  .EO.  6)  GO  TO  9 
K=IV( J) 

S(K)=COS<TH) 

S(K4t )=SIN(TH) 

TH=THtDTH 
1  J=IV(J47) 

IF  (J  .EO.  0)  RETURN 

8  FORNATI^  CIRBND  ERROR', 2110) 

9  URITC(0UTLUN,8>  1,J 
STOP 

END 


SUBROUTINE  CONFL I ( IFLD,CONVAL,N) 

$INCLUDEaRG.DCL)“ 

DIMENSION  CONVALd ),F(3),KV<3),XV<3>,rV(3),XL(2),YL(2) 
J=1 

DO  1  I=1,NUHC 
DO  2  JJ=1,3 
JJJ=J+JJ 
JV=IC(JJJ) 

KA=IV(JV) 

KV( JJ)=KA 
KAIFLD=KA+IFLD 
F( JJ)=S(KAIFLD) 

XV( JJ)=S(KA) *SCALE 

2  rU( JJ)=S(KA+1 )*SCALE 
DO  4  NC=1,N 

NI  =  0 

FC=CONVAL(NC) 

DO  3  JA=1 ,3 
JB=JA+I 

IF  (JB  .EG.  4)  JB-I 

IF  ((F(JA)-FC)*<FC-F(JB))  .LE.  0.)  GO  TO  3 
NI=NI+1 

ALPHAs<FC-FUA)  )/<F<  JB)-F(JA;) 

XL (NI)=(t. -ALPHA )»XV(JA)+ALPHA*XV(JB) 
YL<NI)=(l.-ALPHA)*yV(JA)*ALPHA*YV(JB> 

3  CONTINUE 

IF  <NI  .EO.  0)  60  TO  4 
CALL  LINE<XL,YL,2,1,0,0) 

4  CONTINUE 

I  J=J*ICSIZ 
RETURN 
END 


SUBROUTINE  CQNSRVdFROH) 
♦IMCLUDE<TRG.DCL)' 

SUM=0. 

J=1 

DO  t  l=t,NUMV 
K=IU( J) 

KIFROH=K+IFROH 
SUHaSUN>S(KIFRON) 
t  J=J*IUSIZ 

3  FORHAK'  FIELDM5,"  INTEGRAL',  1F15. 9) 
URITE(0UTLUN,3I  IFRON,SUN 
RETURN 
END 
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SUBROUTINE  COPY( I  TO, IFROM) 
♦INCLUDL(TRG.DCL)' 

J--\ 

DO  I  1=1,NUHV 
K=IV(J) 

KITO=K+ITO 
KIFROH=K+IFRON 
3(KIT0)=S(KIFR0M) 
t  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  CTU(IPOS) 

LOGICAL  NEIGH 
$INCLUDE(TRG.DCL)' 

1  =  1 

DO  1  11=1, NUHV 
NV=IV(I  +  n 
K=iy(I) 

X=S(K) 

Y=S<Ktl) 

KIP0S=K+IP0S 

SIK1F0S)=0. 

DO  2  J=1,NV 
KIPOSJ=KIPOS+J 
2  S(KIP0SJ)=0. 

NVH1=NV-1 
DO  3  JP=1,NVH1 
IJP1=I+JP+1 
NP=IV(IJP1 ) 

KP=IV(NP) 

PX=S(KP) 

Py=S<KP+l ) 

JPP1=JP*1 
DO  3  JO=JPP1,NV 
IJQ1=I+JQ+1 
NQ=1V(IJQ1 ) 

IF  (.NOT.  NEIGH(NP,NQ))  GO  TO  3 
KO=IV(NO) 

QX=S(KO) 

QY=S<KQ+1 ) 

PO»(PX-X)*(OY-Y)-(PY-Y»*(aX-X> 

SIG=1. 

IF  (PO  .LT.  0.)  SIG*-1. 

IF  (IV(JP*I>  .EQ.  6)  00  TO  4 

DX».5*SIG*(PY-Y) 

DY».3*SIG«(PX-X) 

KJP=KIPOS*JP 


iL _ _ _ _ _ _ jt 
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'j(l'.lP0c;)^S'.KlF'US)f.-'’5*DX 
3iKlP0S  +  7)=S(KIP0St,'’)t.75Hir 
S(KJP)-S(KJP)+.25*DX 
S(KJP+7)=S(KJP+7)+.25»DY 

4  IF  (IV(JQ*t)  .EO.  6)  GO  10  S 
DX==-.5*SIG*(QY-Y) 

DY=.5»SIG*(QX-X) 

S(KIP0S)=S(K1P0S)+.75*0X 

S(KIP0S+7)=S(KIP0S+7)+.75*DY 

KJP=KIPOS+JO 

S(KJP)=S(KJP)+.25*DX 

S(KJP+7)=S(KJP+7)+.25*DY 

5  CONTINUE 

DX  =  -SIG»(FY-.5*(QY  +  Y))./3. 
D1=SIG*(PX'.5*(QX+X))/3. 

FX  =  SIG*(aY-.5*(PY»Yn/3. 

EY=-SIG»(QX-.5»(PX*X))/3. 

S(KIP0S)=S(KIP0S)+5.«(DX+EX)/12. 

S(KIP0S+7)=S(K1P0S+7)*5.*(DY+EY)/J 

KJP=KIPOS+JP 

S(KJP)=S(KJP)>DX/6.+5.»EX/12. 

S(KJP+7)=S;KJP*7)+DY/6.+5.*EY/I2. 

KQ=KIPOS+JO 

S(KO)=S(KO)» v.*DX/12.+EX/4. 

S(KQ  +  ,’)=S(K0-7)t5.*0Y/l2.+£Y/A. 

3  CONTINUE 
I=I+IVSIZ 
1  CONTINUE 
RETURN 
END 


SUBROUTINE  CULU(IPOS) 

LOGICAL  NEIGH 
«INCLUDE(TRG.OCL)' 

A1=11./54. 

A2*7./108. 

1  =  1 

DO  1  II=I,NUNV 
HV=IV(I  +  n 
K=IV(I) 

X>S(K) 

Y=S(K+1 ) 

KIPOS=K+IPOS 

S<KIP0S)=0. 

DO  2  J=1,NV 
KIPOSJ=K*IPOS>J 

2  S(K1P0SJ)=0. 

«VH1»HV-1 

DO  3  JP=1,MVH1 
IJP=I+JP 
HP=IV(IJP+1 ) 

KP=IV(«P) 

XPsS(KP) 

YP»S(KP+i ) 

JPP1»JP+1 

DO  3  JO=JPP1,NV 

1JQ«HJ0 

N0=IV(IJQH1 

IF  (.NOT.  NElGH(NP,Na))  GO  TO  3 
KQ>IV(Na> 

XO=S(KO) 

YO>S(Kafl ) 

DELTA*,3*ABS( (XP-X)*(YO-r)-(yP-Y>*(XO-X)> 
KIPOS=K+IPOS 

S(KIPOS)=S(KlPOSHA1  >OELTA 

KIPOSJP«KIPOS*JP 

S(KlPOSJP)>S(KIPOSJP)tA2«OELTA 

KIPOSJO-KIPOStJQ 

S(KIPOSJO)*S<KIPOSja)*A2»DELIA 

3  CONTINUE 
I«I»IVSIZ 

1  CONTINUE 
RETURN 
END 
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SUBROUTINE  DRAUV< IFL , IPUL , IPVL , IP, I t'U, IPV, ASCALE , IPU ) 
tlHCLUDEiTUG.DCL)' 

DIHEMSION  X<5i,V(5) 

£RR=I.E-T 

CALL  UHLUNPtIPU,IPL,lP,ERR) 

CALL  UHUHP<IPU,IPUL,IPU,£RR) 

CALL  UNLU«P(IPU,IPVL,IPV,ERR) 

J=t 

DO  10  1*1,NU«V 
K=IV(J) 

XX=S<K)*SCALE 
YY=S(KH  )*SCAL£ 

KIP=K+lf 

P=S(KIP) 

U=S(KIPU)/P*SCAL£ 

KIPV=KtIPV 

V=S(KIPV)/P»SCALE 

UV=SQRT<U*U+V*V) 

IF  (UV  .LT.  .01»ASCALE)  GO  TO  10 

DX=ASCALE*0*.5 

0Y=ASCALE*V*,5 

Xdl^XX+DX 

X(2)=XX-0X 

Y(1)=YY+DY 

Y<2)=YY-DY 

OX=.4*0X 

DY=,4*0Y 

X(3)*X(1)-DX-DY 

Y(3)*Y(1)-DY*DX 

X(5)  =  X(n-0X*DY 

Y(5)»rU  )-DY-DX 

Y(4)=Y<1) 

X<4)=X<n 

CALL  LINE<X,Y,2,1,0,0) 

CALL  LINE(X<3>,Y<3J,3, 1,0,0) 

10  J=J»IVSIZ 
RETURN 
END 


SUBIVOOTINE  DRCELL 
$INCLUDE(TRG.DCL)“ 

DIMENSION  X(2),Y(:) 

J=) 

DO  1  I^l ,NUNV 
K=IO( J) 

X( I )=S(K)«SCALE 

y(1 )=S(K+1 )*SCALE 

CALL  SYMBOL(X,Y, .05,2,0. 

1  J=J+IVSIZ 
J=1 

DO  2  1=1, NUNC 
JA  =  IC( J  +  1  ) 

JB=IC( J+3) 

JC=IC(J+3) 

KA=IW(JA) 

KB=IVUB) 

t:C=IV(JC) 

XA=S(KA) fSCALE 
YA=S(KAt1 )»SCALE 
XB=S(KB)«SCALE 
YB=S(KB+I )*SCALE 
XC=S(KC)*SCALE 
YC=S(KCt1)*SCALE 
XG=(XA+XB+XC)/3. 
YG=(YA+YB+YC)/3. 

X(1  )=XG 
Yd  )=YG 

X(2)  =  .5«(XA+’XB) 
Y(2)=.5»(YA+YB) 

CALL  LINE(X,Y,2,1 ,0,0) 

X(2)=.5*(XB+XC) 

Y(2)=.5*(YB+YC) 

CALL  LINE(X.Y,2,1 ,0,0) 

X(2)=.5*(XC+XA) 

Y(2)=.5*<YX+YA) 

CALL  LINE(X,Y,2,1 ,0,0) 

2  J=J+ICSIZ 
RETURN 
END 


54 


SUBROUIINE  DPGRID 
$1NCLUDE(TRG.DCL)' 

DIMENSION  X(2),Y(2) 

J=t 

DO  1  I=1,NUHV 
K=IV( J) 

Nv=iv(  j+n 

X(1 )=S(K)»SCALE 

Yd  )=S(CH  )»SCALE 

DO  2  IP=1,NV 

JIP=J+IP 

JP  =  IV( JIP  +  1  ) 

IF  (JP  .6T.  J)  GO  TO  2 
KP-IV! JP) 
X(2)=S(KP)»SCAL£ 
Y(2)=S(KP+n»SCALE 
CALL  L1NE(X,T,2,1 ,0,0) 
2  CONTINUE 
1  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  EUPLOT (!*=•, ZH<^X ,ZMIN) 
«INCLU0E(PA6E.DCL)‘ 

LOGICAL  FAIL 

CHARACTER  ISP, ISTAR, IVERT, IHOR 

DATA  ISP/'  '/, ISTAR/'*'/, IVERT/'iV,IHOR/'_V 

DO  1  J»1,NAPHAXV 

DO  2  I=1,HAPHAXH 

IHAP<I,J)=ISP 

IF  (I  .EQ.  1  .OR.  I  .EQ.  NAPNAXH)  II1AP( I, J)^IVERT 
IF  (J  .EQ.  1  .OR.  J  .EQ.  HAPNAXV)  INAP( I, J)=IHOR 

2  CONTINUE 
I  CONTINUE 

A*2./(HAPHAXH-1 .) 

B=-1.-A 

C=(nAPHAXV-1.)/(Z«IN-ZNAX) 

D=) .-C*ZHAX 
DO  3  1=1, NAPNAXH 
X=A*I+B 

CALL  INTP0L(IP,X,0.,Z,FAIL) 

IF  (FAIL)  60  TO  3 
J=C*Z+D 

IF  (J  .LT.  1  .OR.  J  .6T.  NAPMAXV)  GO  TO  3 
INAP(I,J)=ISTAR 

3  CONTINUE 

4  FORNAT(d ',/,<5X, NAPNAXH  AD) 

URITE(0UTLUN,4)  INAP 

RETURN 

END 


k: 


SUBROUTINE  F LIWV ( 1 FOSfl , I POSB , IPOSC ) 
tIftCLUDEdRB.DCL)' 

J=1 

DO  I  I=1,NU«V 
K=IV( J) 

KIPOSA=K+IPOSA 

KIP0S6=K+IP0fB 

KIPOSC^KUPOS' 

S(KIP0SA)=.5#(S(KIP0SB)+S(KIP0SC)) 
t  J=JtIUSlZ 
RETURN 
END 


SUBROUTINE  FRF ( IPL , IPUL, IPVL , IP, IPU, IPV, IPFX , IPF Y , IPFXL, IPFYL , 
.  DT,F,CF,USX,USY) 

$INCLUDE(TRG.DCL) 

CALL  UNLUNP(2,1PL,IP,1 .E-2) 

CALL  UMLU«P(2,IPUL.IPU,1 .E-2) 

CALL  UNLU«F(2,IPUL,IPV,1 .E-2) 

J=1 

no  I  I=1,NUNV 

IF  (IV(J  +  n  .NE.  6)  GO  TO  1 
K  =  IVU) 

KIP=K*IP 

P=S(KIP) 

KIPU=K+IPU 

U=S(KIPU)/P 

KIPV=K+IPV 

V=S(KIPV)/P 

UV=SQRT(U*U+V»V) 

KIPFX=K+IPFX 

S(KIPFX)=(USX-CFtUfUU) "P 
KIPFY=K+IPFY 

S<KIPFY)=(USY-CF*V«UU)*P 
t  J-JtIVSIZ 

CALL  LUHP(2,IPFX,IPFXL) 

CALL  LU«P(2,IPFr.IPFYL) 

J=1 

DO  2  I=1,NUNV 
K=IV< J) 

An  =  l . 

A12=-DT*F 

AZI'DItF 

A22*1. 

KIPFXL=K+IPFXL 

KIPUL=K+IPUL 

B)=DT*S(KIPFXL)+S(KIPUL) 

KIPFYL=K+IPFYL 

KIPVL=K+IPVL 


56 


B2=PT:»S(KIPFyL)+S(KIPVL) 

D  =  An*A22-Al  1*A21 
S(KIPUL)=<B1*A22-A12»B2)/D 
S(KIPVL)  =  (A1  t*B2-Bl*A2n/D 
2  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  GOCdP) 
$INCLUDE(TRG.DCL) 

1C(NEUC+1 )=IP 

IC(NEyC>2)  =  lFdVSIZ 

lC<NEUC+3)=NEUU-IVSrZ 

IP=IPtIVSIZ 

NEUC  =  NEUCdCSIZ 

RETURN 

END 


SUBROUTINE  GUC(IP) 
lINCLUDETTRe.PCL)' 

IC^NEWC+I )*IP 

IC<NEUC+2>»NEUV 

IC(NEUC+3)=NEUV-IWS12 

NEUV^NEUV^IVSIZ 

NEUC^NEUC^ICSIZ 

RETURN 

END 
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SUBROUTINE  HEXTOP(N) 

»IrtCLUDE(TRG.DCL) 

SCALE=4. 

DO  too  I=1,ICNAX 
too  IC(I)=0 

DO  101  I=1,1VMAX 
101  IV(I)=0 
NEUC=1 
!VSIZ=8 
ICSIZ=4 
IP=t 

NEUV=1+(M+2)*IVSIZ 

IU=N 

DO  1  1=1, M 
DO  2  J=t,IU 
CALL  GUC(IP) 

2  CALL  GDC(IP) 

CALL  GUC(IP) 

IP=IP*IVSIZ 

NEUV=NEUV«IOSIZ 

1  IM=IU+1 
IU=IU-1 
DO  3  1  =  1, N 
DO  4  J=1,IU 
CALL  GDC(IP) 

4  CALL  GUCdP) 

CALL  GDC(IP) 

IP=IP+IVSIZ 

NEUV=NEUVtIVSIZ 

3  IU=IU-1 
NEUV=HEUV-IVSIZ 
NUHC=(N£UC-1 )/ICSIZ 
NUHV=(NEUV-1 )/IVSIZ 
URITE(OUTLUN,S)  NUHC,NUHV 

5  FORNAT('  Allocated  cells  and  vertices' ,2110) 
JC=1 

DO  to  I=t,NUNC 
IVA=IC(JC+1 ) 

IVB=IC(JC+2) 

IVC=IC( JC*3) 

CALL  VLINK<IVA,IVB) 

CALL  VLINK<IVB,IVC) 

CALL  VLINK(1VC,1VA) 

CALL  VLINK(IVB,IVA) 

CALL  VLINK(IVC,IVB) 

CALL  VLINK<IVA,IVC) 

10  JC  =  JCMCSIZ 
IP«1 

DO  20  I=I,NUNV 

IF  (IV(IP+I)  .LT.  6)  GO  TO  21 
20  IP*IP*IVSIZ 


i 

J 
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22  FORHAK'  HO  BOUNDARY  POINTS') 
URlIE(OUrLUN,22) 

STOP 

2t  IB=1P 
IBVF=IB 
IP=0 
NBV=1 

27  IVaB+7)  =  lP 
NV=IV(IB+1> 

DO  23  1=1, NO 
IBI1=lB+I+t 
J=IV(IBI1 ) 

IF  (J  .£0.  IP)  GO  10  23 
IF  (J  .EQ.  IBVF)  GO  TO  25 
IF  (I0(J*1)  .LT,  6)  GO  TO  24 

23  CONTINUE 

26  FORNAK'  BOUNDARY  ERROR') 

URITE(0UTLUN,26) 

CALL  OCDUHP 
STOP 

24  NBV=NBU41 
IP«1B 
IB=J 

GO  TO  27 

30  FORHAK'  NUMBER  OF  BOUNDARY  VERTICES',!  10) 

25  URITC(0UrLUN,30)  NBV 
IBVF=IB 

RETURN 

END 


SUBROUTINE  INITF (IP, IPU , IPV ) 
«INCLU0E(TRG.DCL)^ 

.1=1 

DO  1  I=I,NUNV 
K=IV( J) 

KIP*KtIP 
S(KIP)  =  1 . 

KIPU=K+IPU 

S(KIPU)=0. 

KIPV»K+IPV 

S(KIPV)=0. 

I  J=J+IVS1Z 
RETURN 
END 


4 
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SUBROUTINE  INTPOK IFLO,X, r,2.FAIL> 

♦INCLUDE(TRG.DCL)‘ 

LOGICAL  FAIL 
FAIL=.TRUE. 

ICELL=1 

DO  to  I=1,NUHC 
IP=IC(ICELL  +  1 ) 

IQ=IC(ICELL+2) 

IR=^IC(ICELL  +  3) 

JP  -IVdP) 

JQ=IV(IQ) 

JR=IV(1R) 

XP=S( JP)-X 

YP=S(JP+1 )-Y 

XQ=S(JQ)-X 

YQ=S(JQ+1 )-Y 

XR=S(JR)-X 

YR=S( JR+1 )-Y 

PQ=XP*YO-YP»XQ 

QR»XQ*YR-YQ»XR 

RP=XR*YP-YR*XP 

IF  (Pa*OR  .LT.  0.)  GB  TO  1 

IF  <QR*RP  .L7.  0.)  GO  TO  1 

IF  (PQ*RP  .LT.  0.)  GO  TO  1 

JPIFLD=JP+IFLD 

ZP=S(JPIFLD) 

JOIFLD=JOdFLD 

ZO»S<JQIFLD) 

JRIFLDsJR+IFLD 

ZR=S(JRIFLD) 

0=XP* ( YQ- YR ) - YP M  XO -XR ) ♦Xa»YR- YO^XR 

Z=<XP*(YQ*ZR-ZO»YR)-YP*(XO*2R-ZO=»XR)+ZP*(XO*YR-YO*XR))/D 
FAIL=. FALSE. 

RETUFL 

1  ICELL=ICELL+ICSIZ 
to  CONTINUE 
RETURN 
END 


I 


FUNCTION  INTRI(X,Y,ICELL) 
♦INCLUOE(TRG.DCL)‘ 

LOGICAL  INTRI 
IP=IC(ICELL*1) 
IQ=IC(1CEI.L*2) 
IR=IC(ICELL+3) 

JP=IV<IP) 

JQMV(IQ) 

JR=IV<JR) 

XP=S(JP)-X 

YP=S(JP+1 >-Y 

XO=S(JQ)-X 

YR«S(JQ+1)-Y 

XR=S( JR)-X 

YR=S<JR+1 )-Y 

PQ»XP*YO-YP»Xa 

QR=XQ«YR-YO«XR 

RP=XR*YP-YR*XP 

IF  (PO*QR  .LT.  0.)  GO  TO  1 

IF  (QR*RP  .LT.  0.)  GO  TO  1 

IF  (PQ*RP  .LT.  0.)  GO  TO  1 

IHTRI=.TRUE. 

RETURN 

1  INTRI*. FALSE. 

RETURN 

END 
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SUBROUTINE  LUMP ( IPOS, IFROM, ITO) 
«INCLUDE(TRG.DCL)‘ 

J=1 

DO  1  I=1,NUHU 
K=IV(J) 

NV=IV(J+1) 

KITO=K*ITO 

KIPOS=K+IPOS 

KIFROH=K>IFRON 

S(KITO)=S(KIPOS)*S(KIFROM) 

DO  2  IN=I,NV 
JINl=J+IN+t 
JP^UMJINI ) 

KP=IV( JP) 

KIPIN=KIPOS*IN 

KPIFR=KP+IFRON 

2  S(KITO)=S(KITO)+S(KIPIN)*S(KPIFR) 
1  J=iJ*IVSIZ 
RETURN 
END 


SUBROUTINE  HAXNIN(IFLO) 
(INCLUOE(TRG.DCL)* 

XNAX«S<IFLO*1 ) 

XNIN^XNAX 

J=1 

DO  1  i=i,Nunv 
K=IV( J) 

KIFLD=K+IFLD 

X=S<KIFLD) 

IF  (X  .GT.  XHAX)  XNAX=X 
IF  <X  .LT.  XNIN)  XNIN=X 

1  J=J+IVSIZ 

2  FORNAK'  nAX,niN'',2FI0.6) 
URITE<0UTLUN,2)XHAX,XHIN 
RETURN 

END 


FUNCTION  NEIGHdP.IQ) 
$INCLUDE<TRG.DCL)' 

LOGICAL  NEIGH 
NV=IV(IP+1 ) 

DO  1  1=1, NV 
IPI1=IPtI+1 

IF  (IV(IPII)  .EQ.  IQ)  GO  TO  2 

1  CONTINUE 
NEIGH=. FALSE. 

RETURN 

2  NEI6H=.TRUE. 

RETURN 

END 


SUBROUTINE  PRCON< IFLP,XNIN,XMAX) 
LOGICAL  FAIL 

DINENSION  ICON(t3),LIN(nAPHAXH) 
DATA  icon/'-' ,'Q- ,'2' ,'2' ,'A' 
BX=-(NAPHAXH«1 . I/HAPNAXH 
BT=(nAPHAXV*1.>/HAPNAXU 
URITE(OUTLUN,5) 

DO  1  J=t,HAPNAXV 

Y=-J/(HAP«AXV/2.)*BY 

DO  2  I=t,HAPHAXH 

X=I/(«AP«AXH/2.)+BX 

R=X*X+Y*Y 

K=13 

IF  (R  .GT.  1.)  GO  TO  3 
CALL  INTPOL(IFLD,X,Y,Z,FAIL) 

IF  (FAIL)  GO  TO  3 
K=10.*(2-XHIN)/(XHAX-X«IN) 

IF  (K  .LT.  0)  K=--1 
IF  (K  .GT.  10)  K=10 
K=K*2 

3  LIN(I)=ICON(K) 

2  CONTINUE 

4  F0RHAT(3X,NAPNAXH  At) 
URITE(0UTLUN,4)  LIN 

1  CONTINUE 

5  FORHAT(d') 

RETURN 

END 
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SUBROUTINE  PRFLDSdt  ,12) 
tlNCLUDE(TRG.OCL)' 

DINENSION  VdO) 

URITE(0UTLUN,4)  11,12 
4  FORHftK'  FIELDS', M,'  THROUGH',14) 
URITE(0UTLUN,4)  It, 12 
.J=1 

DO  1  I<>1,NUHV 
K=IV(J) 

NN=1 

DO  2  N=I1,I2 

KN=KtN 

V(NN)=S(KN) 

NN=NN+1 

IF  (NN  .LE.  10)  GO  TO  2 
URITE(0UTLUN,3)  I,V 
3  FORNRTI'  ',I8,10F12.A) 

NN=1 

2  CONTINUE 
NN=NN-1 

IF  (NN  .GT.  0)  URITE<0UTLUN,3)  I,(U(N),N>1,NN) 
1  J*J+IVSI2 
RETURN 
END 


SUBROUTINE  QPRCON ( IFLD, XMI N , XHAX ) 
CHARACTER  IC0N(13) 

«INCLUDE(TRG.DCL)' 

«INCLUDE(PAGE.DCL)' 

DATA  IC0H/'-','0','1','2','3','4','5 
DO  2  J=1,NAPNAXV 
DO  2  I-1,NAPHAXH 
2  INAP<I,J)-IC0N(13) 

IVERT=1 

DO  I  II=I, NUNV 
KK::IO(IVERT) 

X=S(KK) 

Y=S(KK  +  1 ) 

KKIFLD=KK+IFLD 

Z=S(KKIFLO) 

I*  (HAPNAXH/2.)*X*HAPHAXH/2.*I. 

IF  (I  .LT.  1)  1=1 

IF  (I  .GT.  NAPNAXH)  I>HAPHAXH 

J*-(HAPHAXV/2. )*Y*HAPHAXV/2.+1 . 

IF  (J  .LT.  1)  J=1 

IF  <J  .GT.  NAPHAXV)  J=NAPHAXV 

K=tO.«(Z-XNIN)/(XHAX-XNIN) 

IF  (K  .LT.  0)  K=-I 
IF  <K  .OT.  10)  K>10 


+•' , 


V 


■i 


64 


K=K*2 

IHAP(I,J)=ICON(K) 

I  IVERT=IVERT+IWSIZ 
3  F0RHAT(M%/,(3X,HAPHAXH  AD) 
URITE(0UTLUN,3)  IHAP 
RETURN 
END 


SUBROUTINE  RIGBNDdPUL,  IPVL,  IPU,  IPV) 
$INCLUDE(TRG.DCL)‘ 

CALL  UNLUNP(2,IPUL,IPU,1.E-3) 

CALL  UNLUNP(2,IPVL,IPV,l.E-5) 

J=IBVF 

DO  1  1=1, NBV 
K=IV( J) 

KIPU=K+IPU 

S(KIPU)=0. 

KIPV=K+IPV 

S<KIPV)=0. 

1  J=IV(J*7) 

CALL  LUHP(2,1PU,IPUL) 

CALL  LUNP<2,IPU,IPVL) 

RETURN 

END 


SUBROUTINE  SBROKOHEGA,  IPOSU) 
«INCLUDE<TRG.DCL)' 

1  =  1 

DO  1  11=1, NUNV 
K=IV(I) 

X=S(K) 

Y=S(K+D 
I=I+IVSIZ 
KIPOSU=KtIPOSU 
S<KIPOSU)>-Y*ONEGA 
1  S(KIPOSU>1)>X*ONEGA 
RETURN 
END 
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SUBROUTINE  STEP(IPLN,IPULN,IFVLN,IPL0.1fUL0,IPVL0.IPL,IPUL,IPVL, 
.  IP,lPU,IPV,IPUU,IPuO,IPVV,lPFX,iPFY,lPFXL,iPFYL,lPU,DT, 

.  F,CF,WSX,WSY) 

«1NCLUDE(TRG.0CL) 

CALL  0NLUHP(2,IPL,IP,t.E-2) 

CALL  UHLUHP<2,IPUL,IPU,1.E-2) 

CALL  UNLUHP(2,IPyL,lPV,l .E-2) 

J=1 

DO  1  I=1,NUNV 
K=IV(J) 

KIP=K+IP 

KIPU=K+IPU 

KIPV=K*1PV 

KIPUU=K+IPUU 

KIPUV=K+IPUV 

KIPVg=K+IPVV 

P=S<KIP) 

U=S<KIPU)/P 

V=S<KIPV)/P 

S(KIPUU)=P*(U*U+.5=»P) 

S<KIPUV)=P*U*V 

S<KIPVg)=P*<V*V+.5*P) 

1  J»J+1VSIZ 

CALL  TRNSPT<IPLN,IPLO,IPU,IPV,IPU,DT) 

CALL  TRNSPT(IPULN,IPULO,IPUU,IPUV,IPU,DT) 

CALL  TRHSPT(IPVLN,IPMLO,IPUV,IPVM,IP«,DT) 

CALL  FRF(IPLH,IPULH,1PVLH,IP,1PU,1PV,IPFX,1PFY,IPFXL,IPFYL,DT, 

.  F,CF,USX,USY) 

CALL  RIGBNO<IPULN,IPVLN,IPU,IPV> 

RETURN 

END 


SUBROUTINE  SUGRID(MITER) 
♦INCLUDE<TRG.DCL)* 

RELPA=1.388 
DO  9  NI  =  1, NITER 
J=1 

CHNAX=0. 

DO  I  I=t,NUHO 
NV=IV(J+1 ) 

IF  (NV  .LT.  6)  GO  TO  1 
X=0. 

Y*0. 

DO  2  K=1,& 

JK1=J+1+K 
JN>IV( JK1 ) 

JS=IV(JN) 

X=X+S(JS) 

2  Y=Y+S(JS*1) 

JS=IV( J) 

XO-S< JS) 

YO=S(JS+l ) 

X=X/6. 

Y»Y/6. 

CH*SORT( ( X-XO) ►•2* ( Y-Y0)**2) 

IF  (CH  .GT.  CHHAX)  CHHAX^CH 
S<JS)»RELPA*X*(1 .-RELPA)»S(JS) 
S(JS+1)=RELPA*Y+(1-RELPA)*S(JS+1 ) 

1  J»J+IVSI2 

8  FORHAK"  HAXINUn  DISPLACEMENT', IFI2. 8) 
URITE(0UTLUN,8}  CHHAX 

9  CONTINUE 
RETURN 
END 


rii 


3UBR0U1 INE  TRNSF'l  ( IN!- U  .  lOLf,  I(  X,  IF  t ,  IPU.  DT  ) 
SINCLUHEdRG.DCL) 

J=I 

DO  t  I=1,NU«V 
K=IV(J) 

NV^IV( J+l ) 

K1NEU=K+INEU 

KIFX=K+IFX 

KIPU=K+IRU 

KIFY=K+IFY 

S(KINE«)=S<KIFX)*S(KIPU)+S(KIFY)*SlKIPU+7) 

DO  2  IN=1 ,NV 
JIN=J+IN 
JP=IV( JIN+l ) 

KP=IV( JP) 

KPrFX=KP+IFX 

KIPUIN=KIPU+IN 

KPIFY=KP+IFY 

2  S(KINEU)=S(K1NEU)+S(KPIFX)*S<KIPU1N)+S(KPIFY) xSCKIPUINt?) 
KIOLD=K+IOLD 

S(K1NEU)=:S(KI0LD)-DT*S(KINEU) 

1  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  TSTCTU(IPOS) 

*INCLUDE(TRe.DCL)' 

URITE<OUTLUN,<) 

A  FORNRK"  TEST  T-UEIGHTS  V8X , '  I '  ,8X,  NV' ,  1 8X , 'TEST  X',HX,'TEST  Y  ) 
1  =  1 

DO  I  11=1, NUNV 
NV=IVn  +  1  ) 

K=IV(I) 

KIPOS=KMPOS 

X=S(KIPOS) 

Y=S(KIPOS+7) 

DO  2  J=1 ,NV 

KIPOSJ=KIPOS+J 

X=X+S(KIPOSJ> 

Y=Y*S(KIP0SJ+7) 

2  CONTINUE 

3  F0R«AT<1X,2II0,2F20.10) 

URITE(OUTLUN,31  I,ND,X,Y 

I  I=I>IVSIZ 
RETURN 
END 


hf 


SUFi|-:OUi  INI;  UNLUrtPl  IPOS.IFI<Of1,ITO,LRR) 
»I«CLUIi£(  IRG.DCL)  * 

RtLPA^I  .3 
N1TER=25 

DO  1 1  Nir^l  .NITER 
CH«AX=0. 

J=1 

DO  10  I  =  1,NUI1V 
K=IV( J) 

KITO=K+ITO 
SOLD=S(KITO) 
t;lFROh  =  K+IFRON 
S(KITO)=S(KIFROn) 

NV=IV(  JH  ) 

DO  12  IN=1 .NO 
JIN^J+IN 
JP=^IV(  JIN  +  1  ) 

KP=IV( JP) 

KIPSIN=K+IPOS+IN 

KPIIO=KP+ITO 

12  S(KITO)=S<KITO)-S(KIPSIN)*S(KF'ITO' 
KIPOS=K+IPOS 

3(KIT0)=S(KIT0)/S(KIP0S) 
CH=ABS(S(KITO)-SOLD) 
S(KITO)=RELPA*S(KITO)+( 1 .-RELPA )»SOLD 
IF  (CHHAX  .LT.  CH)  CHHAX=CH 

10  J=J+IVSIZ 

IF  (CHHAX  .LT.  ERR)  GO  TO  H 

11  CONTINUE 

13  FORHATl  HAXinUH  CHANGE  IN  UNLUMP  ,IF12. 8) 
URITE(0UTLUN.13)  CHHAX 

M  CONTINUE 
RETURN 
END 


fi9 


SUBKOUTINL  VCDUflF' 

IINCLUCEITRG.DCL) 

URITE(OUTLUN, 1  ) 

1  FORMATC  CELL  LIST') 

J=l 

DO  2  1=1, NUNC 
JL=J+ICSIZ-1 

URITE<0UTLUN,3)  I , J , ( IC ( K ) ,K=J, JL) 

2  J=J+1CSIZ 

3  F0RNAT(1X,6I10) 

J=1 

4  FORNAK'VERTEX  LIST') 
URITE<0UTLUN,4) 

DO  5  I=1,NUNV 
JL=J+IVSIZ 

URITE(0UTLUN,6)  I ,  J, ( IV(K ) ,K=J, JL ) 
6  F0RMAT(<1X, 10110)) 

5  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  VLINKdVA,  IVB) 

♦  INCLUDEdRG.DCD* 

NV  =  IV(IVA  +  1 ) 

IF  (NV  .EQ.  0)  GO  TO  1 
DO  2  1=1, NV 
IVAII=IVA*I+1 

IF  (IV(IVAII)  .EG.  IVB)  RETURN 
2  CONTINUE 
1  IVdVA  +  l  )=NV  +  1 
IVANV=IVA+NV 
IVdVANV+2)=IVB 
RETURN 
END 


SUBROUTINE  VST0d0RI6,NUDS) 
lINCLUDEdRG.DCD* 

K=IORIG 

J=l 

DO  1  I=I,NUHV 
IV(J)=K 
K=K+NUDS 
1  J=J+IVSIZ 
RETURN 
END 


END 

DATE 
FILMED 


DTIC 


